home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / byte0987.arc / KAREX.ARC / KAREX2.BAS < prev    next >
Encoding:
BASIC Source File  |  1980-01-01  |  4.0 KB  |  135 lines
  1. 100 '-------------------------------------------------------
  2. 101 '
  3. 102 '   KAREX2.BAS is a Microsoft BASIC Release 5 program
  4. 103 '   that solves EXAMPLE 2 of the article
  5. 104 '
  6. 105 '              KARMARKAR'S ALGORITHM
  7. 106 '
  8. 107 '     by Andrew M. Rockett and John C. Stevenson
  9. 108 '
  10. 109 '    This program was written by Andrew M. Rockett
  11. 110 '
  12. 111 '-------------------------------------------------------
  13. 200 '
  14. 202 ' N is number of unknowns and K is the number
  15.       of equations
  16. 204 '
  17. 206 N = 8 : K = 4
  18. 208 '
  19. 210 K1 = K + 1 : K2 = 2*K1
  20. 212 DIM A0(N), XOLD(N), XNEW(N), CC(N), CP(N),
  21.     A(K,N), B(K1,N), B1(K1,K2), B2(N,K1), B3(N,N)
  22. 214 FOR C = 1 TO N : A0(C) = 1/N : XNEW(C) = A0(C) :
  23.     NEXT C
  24. 216 '
  25. 218 ' T is the tolerance
  26. 220 '
  27. 222 T = .001
  28. 224 '
  29. 226 ' ALPHA is usually set equal to 1/4 ...
  30. 228 '
  31. 230 ALPHA = .25
  32. 232 '
  33. 234 ITERATION = 0
  34. 236 '
  35. 238 ' Data for constraint matrix A
  36. 240 '
  37. 242 DATA 1,  0,  1,  0,  0,  0,  1, -3
  38. 244 DATA 1,  0,  0, -1,  0,  0,  2, -2
  39. 246 DATA 0,  1,  0,  0,  1,  0,  3, -5
  40. 248 DATA 0,  1,  0,  0,  0, -1,  4, -4
  41. 250 '
  42. 252 FOR R = 1 TO K : FOR C = 1 TO N : READ A(R,C) :
  43.     NEXT C : NEXT R
  44. 254 '
  45. 256 ' Data for objective function CC
  46. 258 '
  47. 260 DATA 0,  0,  0,  0,  0,  0,  1,  0
  48. 262 '
  49. 264 FOR C = 1 TO N : READ CC(C) : NEXT C
  50. 266 '
  51. 268 V = 0 : FOR C=1 TO N : V = V + CC(C)*A0(C) : 
  52.     NEXT C : VNEW = V
  53. 270 '
  54. 272 ' Main iteration process is the same as in 
  55.       KAREX1.BAS ...
  56. 274 '
  57. 300 WHILE VNEW/V > T
  58. 301 PRINT USING "###";ITERATION;:
  59.     FOR C=1 TO N:PRINT USING "###.####";XNEW(C); :
  60.     NEXT C :
  61.  PRINT USING "####.#######";VNEW/V
  62. 302 ITERATION = ITERATION + 1
  63. 303 FOR C = 1 TO N : XOLD(C) = XNEW(C) : NEXT C
  64. 304 FOR R=1 TO K:FOR C=1 TO N:B(R,C)=A(R,C)*XOLD(C):
  65.     NEXT C:NEXT R
  66. 305 FOR C=1 TO N:B(K1,C)=1:NEXT C
  67. 306 FOR R=1 TO K1 : FOR C=1 TO K2 : B1(R,C)=0 :
  68.     NEXT C : NEXT R
  69. 307 FOR R=1 TO N  : FOR C=1 TO K1 : B2(R,C)=0 :
  70.     NEXT C : NEXT R
  71. 308 FOR R=1 TO N  : FOR C=1 TO N  : B3(R,C)=0 : 
  72.     NEXT C : NEXT R
  73. 309 FOR C=1 TO N  : CP(C) = 0 : NEXT C
  74. 310 FOR R=1 TO K1:FOR C=1 TO K1:
  75.      FOR I=1 TO N:B1(R,C)=B1(R,C)+B(R,I)*B(C,I):
  76.       NEXT I:
  77.     NEXT C:NEXT R
  78. 311 FOR I = 1 TO K1 : B1(I,I+K1)=1 : NEXT I
  79. 312 FOR R = 1 TO K1
  80. 313  IF B1(R,R) <> 0 THEN 318
  81. 314    I = R + 1
  82. 315    IF I > K1 THEN PRINT "Error! BBT is SINGULAR!" : 
  83.        GOTO 407
  84. 316    IF B1(I,R) = 0 THEN I = I+1 : GOTO 315
  85. 317    FOR C = 1 TO K2 : SWAP B1(R,C),B1(I,C) : NEXT C
  86. 318  FOR I = R+1 TO K1:Z = B1(I,R)/B1(R,R):
  87.       FOR C=1 TO K2:B1(I,C)=B1(I,C)-Z*B1(R,C):NEXT C:
  88.      NEXT I
  89. 319 NEXT R
  90. 320 FOR R=K1 TO 2 STEP -1:FOR I = R-1 TO 1 STEP -1:Z =
  91.      B1(I,R)/B1(R,R):
  92.      FOR C=R TO K2:B1(I,C)=B1(I,C)-Z*B1(R,C):NEXT C:
  93.     NEXT I:NEXT R
  94. 321 FOR R=1 TO K1:Z = B1(R,R):
  95.      FOR C=1 TO K2:B1(R,C)=B1(R,C)/Z:NEXT C:
  96.     NEXT R
  97. 322 FOR R=1 TO N:FOR C=1 TO K1:
  98.      FOR J=1 TO K1:B2(R,C)=B2(R,C)+B(J,R)*B1(J,C+K1):
  99.       NEXT J:
  100.     NEXT C:NEXT R
  101. 323 FOR R=1 TO N:FOR C=1 TO N:
  102.      FOR J=1 TO K1:B3(R,C)=B3(R,C)+B2(R,J)*B(J,C):
  103.       NEXT J:
  104.     NEXT C:NEXT R
  105. 324 FOR R = 1 TO N : B3(R,R) = B3(R,R) - 1 : NEXT R
  106. 325 FOR R=1 TO N:FOR C=1 TO N:B3(R,C)=-1*B3(R,C):
  107.     NEXT C:NEXT R
  108. 326 FOR R=1 TO N:FOR C=1 TO N:B3(R,C)=B3(R,C)*XOLD(C):
  109.     NEXT C:NEXT R
  110. 327 FOR R=1 TO N:FOR C=1 TO N:CP(R)=CP(R)+B3(R,C)*CC(C):
  111.     NEXT C:NEXT R
  112. 328 AA=0:FOR C=1 TO N : AA = AA + CP(C)*CP(C) : NEXT C
  113. 329 AA = SQR(AA) : FOR C=1 TO N : CP(C) = CP(C)/AA : 
  114.     NEXT C
  115. 330 AA = SQR(N*(N-1))/ALPHA
  116. 331 FOR C=1 TO N : XNEW(C) = (A0(C) - CP(C)/AA)*XOLD(C) : 
  117.     NEXT C
  118. 332 AA=0:FOR C=1 TO N : AA = AA + XNEW(C) : NEXT C
  119. 333 FOR C=1 TO N : XNEW(C) = XNEW(C)/AA : NEXT C
  120. 334 VNEW=0:FOR C=1 TO N:VNEW=VNEW+CC(C)*XNEW(C):NEXT C
  121. 335 WEND
  122. 336 '
  123. 400 PRINT:PRINT"Tolerance reached: Vnew/Vinitial = ";
  124.     VNEW/V:PRINT
  125. 401 PRINT USING "###"; ITERATION; :
  126.     FOR C=1 TO N:PRINT USING "###.####";XNEW(C); : NEXT C :
  127.     PRINT USING "####.#######";VNEW/V
  128. 402 '
  129. 403 ' Project solution from simplex back to orthant ...
  130. 404 '
  131. 405 PRINT:FOR C=1 TO N-2:PRINT XNEW(C)/XNEW(N), :
  132.     NEXT C:PRINT
  133. 406 '
  134. 407 END
  135.